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PREFACE 


This unit was written to give the student of statistical mechanics an 
opportunity to bring the power of the computer to bear on this fascinating 
subject. The computer can be used as a very effective tool to probe the 
concepts involved. 


It is assumed that this unit is being used as a supplementary text in a 
course in statistical mechanics. The only requirements besides being en- 
rolled in or having taken such a course is some familarity with BASIC 
language programming and access to a computer. 


Many of the exercises involve writing or modifying computer programs. 
This is intentional and follows from the strategy that one learns more and 
becomes more fundamentally involved by writing programs than by 
merely interacting with programs already prepared. The required level of 
skill in programming is not high. Students not familar with BASIC should 
be able to do fine after a couple of hours of instruction or study on their 
own. 


Students and instructors as well should avoid falling into the trap of 
thinking that here is a right or best method to solve a specific problem. 
Certainly some programs may be more elegant or efficient than others. 
However, to become obsessed with intricate solutions or flashy details of 
programming is to draw attention away from the physics. Remember that 
the purpose of this unit is not to teach about computers, but to use 
computers to learn about physics. 


This unit is not meant to be a text, and that no effort has been made to 
provide a complete treatment of statistical mechanics. Excellent texts are 
available and should be used for completeness. This unit is designed to 
reinforce such a text and follows closely the approach used by Reif? and 
also draws heavily upon the unique description of statistical mechanics by 
Sherwin‘. 
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CHAPTER ONE: INTRODUCTION 


MACROSTATES AND MICROSTATES 


We will begin this unit with a brief definition of macrostates and micro- 
states since this is at the heart of statistical mechanics. The prefixes macro 
and micro are defined with respect to atomic dimensions. A microscopic 
object has dimensions of the order of magnitude of the atomic diameter 
(about 10° cm). On the other hand, a macroscopic object has dimensions 
very much larger than the atomic size. Macroscopic objects are ‘‘people- 
sized’’ and constitute the everyday world of our experience. 


By macrostate, we mean a complete description of a macroscopic system 
using macroscopic parameters. For example, we might describe a container 
of gas in terms of its volume, pressure, and temperature. These are large 
scale measurements which are readily available to us. A microstate is, as we 
would expect, a description of a microscopic system using microscopic 
parameters. For example, we could specify the mass, position, and velocity 
of each of the particles in the container of gas at some instant. It is 
obvious that this very large quantity of microscopic information is gen- 
erally not available to us. However, it is also clear that in some fashion the 
macrostate is functionally dependent upon the microstate. 


Our fundamental problem in statistical mechanics is to find a way to 
understand the behavior of macroscopic systems containing very large 
numbers of particles (of the order of magnitude of Avogadro’s Number, or 
about 107°) by reasoning from known first principles at the atomic level. 
It is easy to see that this is hopeless if approached from the point of view 
of classical mechanics. Given one particle moving under the influence of 
known forces, it is a routine matter to determine its motion and momen- 
tum as functions of time. With a computer, we can do the same thing for 
several particles. However, with 102% particles it is clear that the task 
cannot be accomplished with the techniques of classical mechanics, even 
with the largest computer. Another approach is required. 


The necessary bridge between the microstate and corresponding macro- 
state is provided by statistical mechanics. This is a very well understood 
and complete branch of physics that has been developed in the last 100 
years. With statistical mechanics we can, in fact, predict descriptions of 
macrosystems from microscopic information. It is thus interesting to note 
that solutions are available for the two extremes — one particle, or very 
many particles — but that both classical mechanics and statistical me- 
chanics fail in the region between. 


HISTORICAL PERSPECTIVE 


Quite often, the terms thermodynamics and statistical thermodynamics are 
carelessly substituted for statistical mechanics, causing confusion in the 
minds of many students. From a historical point of view, the differences 
between these terms are interesting and provide valuable insight. 


Thermodynamics is generally understood to be defined by the following 
four laws: 


1. Zeroth Law 


Two systems, each of which is in equilibrium with a third system, 
are in equilibrium with each other. 
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2. First Law 
AE=W+Q, 


where AE is the change in mean energy of a system, W is the work 
done on the system, and Q is the heat transferred to the system. 


3. Second Law 
ds = and AS>=0, 


where S is the entropy of a system, and T is the absolute tempera- 
ture at which an amount of heat dQ is transferred to the system. 


4, Third Law 
AsT > 0, 3 So j 


where Sq is a constant independent of the system. 


These laws which define thermodynamics are completely macroscopic in 
nature. They were discovered and applied long before a comprehensive 
theory of atoms and molecules was available. Note that the laws make no 
reference to the detailed structure of the system, nor ts there any element 
of probability present. Functionally, however, thermodynamics is a useful 
tool even though it cannot bridge the gap between the microstate and 
macrostate. 


If we add the notion of probability to thermodynamics, we define statis- 
tical thermodynamics. Probability is introduced through the statistical 
relation 


Pa eS/k , 


where P is the probability of observing a macrostate with entropy S, and k 
is the Boltzmann constant. 


Finally, if we connect the entropy S to the number of microstates 
accessible to the system, we have statistical mechanics. Statistical me- 
chanics is therefore the broadest subject, with statistical thermodynamics 
and thermodynamics being progressively more narrowly defined. The key 
to statistical mechanics is the relation 


S = k In(W) 


which is credited to the German physicist Ludwig Boltzmann (1844- 
1906). The significance of this relation is underlined by the fact that it is 
carved into the Boltzmann memorial in the Central Cemetery in Vienna. A 
more modern version is 


S = k In({2) 


where k is the Boltzmann constant and {2 is the number of microstates 
accessible to the system. Josiah Williard Gibbs (1839-1903), the first 
noteworthy American theoretical physicist, used this modern version and 
the concept of statistical ensembles to construct a very general and 
beautiful framework for statistical mechanics that remains valid in the 
most modern applications. 





THE ROLE OF THE COMPUTER 


It is important to establish the role of the computer with respect to this 
unit on statistical mechanics. Most of the essential ideas of statistical 
mechanics can be demonstrated by examining the behavior of systems 
containing from 20 to 100 elements. The computer is required to carry 
out the lengthy calculations which usually would preclude this type of 
investigation. Consequently, the computer is a tool which gives us the 
leverage we need to get at the ideas of statistical mechanics. Using the 
computer as a device to simulate the random processes involved, we can 
uncover those ideas which have fundamental importance. /t is precisely 
those ideas which are most often slighted in favor of a cookbook-formula 
approach to statistical mechanics. 
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CHAPTER TWO: PROBABILITY, 
STATISTICAL MEASUREMENT, AND 
RANDOM PROCESSES 


The concepts of probability and random processes are fundamental to any 
modern treatment of statistical mechanics. Systems reacting to the influ- 
ence of blind chance behave in predictable ways. Consequently, we must 
carefully consider the topics of probability and statistical measurements If 
we hope to understand the behavior of statistical systems. 


LAWS OF PROBABILITY 


We will define the laws of probability in terms of states. By state we mean 
a complete description of a system or event. The fundamental law of 
probability defined in terms of states is 


Number of states which imply A 


P(A) = 
nN Total number of possible states 


(1) 


P(A) is read as ‘‘The probability of A.’ If all the states imply A, then P(A) 
= 1. If none of the states imply A, then P(A) = 0. Thus, al! probabilities 
fall in the range 0 to 1 inclusive. 


Two examples will illustrate the fundamental law: 


1. If we deal a card from a well-shuffled bridge deck of cards, what is 
the probability of getting a King? 


P(K) = Number of states corresponding toa King 4 
Total number of possible states 52 - 


2. If a pair of dice is thrown, what is the probability that the faces 
uppermost will add to 7? In this case we can describe a state with a 
number pair. The first number is the uppermost face on die 1, the 
second number refers to die 2. Thus, (3,1) indicates a3 on die 1, and 
a 1 on die 2. Note that (3,1) and (1,3) are different states. It is easy 
to enumerate the states which correspond to a sum of seven. They 
are (1,6), (6,1), (2,5), (5,2), (3,4), and (4,3) for a total of 6 states. 
Since to each of the six faces of die 1 we could have six faces of die 
2, there are a total of 36 outcome states possible. Thus 


6 
f7)=— 
P(Sum of 7) 36 


We can use (1) to define other laws of probability. Suppose we want to 
Know the probability of e/ther event A or event B. The law which gives us 
this probability is 


P(A or B) = P(A) + P(B) . {2) 


Returning to the deck of cards, we can easily compute the probability of 
dealing either a heart or a club: 


P(Heart or Club) = P(Heart) + P(Club) 
= 13/52 + 13/52 = 26/52. 
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To follow the idea expressed by (2) to an important conclusion, suppose 
that A, B, C, and D are the only possible outcome states for a system. It 
follows that 


P(A or B or C or D) = P(A) + P(B) + P(C) + P(D) = 1 


This is equivalent to the statement that something must happen. We 
generalize this concept to 


P(One of the possible outcome states occurring) = 1 (3) 


Another type of situation we must be able to handle is that of an event A 
followed by an event B. If the two events are independent (the occurrence 
of the second event is not related to the occurrence of the first) we have 


P(A and B) = P(A)P(B) . (4) 


However, if the two events are dependent (the occurrence of the second 
event is related to the occurrence of the first) the desired relation is 


P(A and B) = P(A}P(BIA) . (5) 
P(BIA) is read as ‘’The probability of B given the occurrence of A.” 


some examples should clarify the difference between dependent and 
independent events. 


1. If a deck of cards is shuffled, a card is dealt, and then returned to 
the deck which is again shuffled. A second card is dealt. What is the 
probability of getting a spade followed by a King? Clearly the two 
events are independent since the occurrence of the first event cannot 
affect the second. Therefore we use (4), 


P(Spade and King) = P(Spade)P(King) 
= (13/52) (4/52) = 1/52 


2. Now suppose that two cards are dealt in succession. Now the events 
are dependent since the occurrence of the first event bears on the 
second, and we should use (5). What is the probability of getting the 
Ace of hearts, then any other heart? 


PlAnearts and Heart) = PlAnearts)P(HeartlAy garts) 
= (1/52)(12/51) = 12/2652 


A final idea concerning probability which is needed involves combinations. 
For example, suppose we want to know the probability of getting two 
heads and three tails if five coins are flipped at once. Any of the following 
states correspond to our desired event: 
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Therefore, we want 
P(State 1 or State 2or * * * or State 10). 


However, it can be seen that 


P(State 1) = P(State 2}= «+ « * =P(State 10). 


Consequently, the probability of two heads and three tails when five coins 
are tossed is given by 


P(2H and 3T) = 10 P0(HRH&T&T&T) = 10 (1/2)5. 


In this case it was easy to enumerate all the possible states. In general, 
however, this is not true and we must be able to compute the number of 
states. If N. is the number of states corresponding to N items taken n ata 
time, then 


| 
y =--__N! 


r onl(N- on)!” (6) 


If N=5, andn=2 


N-=sar = ayy 7 10, 


which agrees with our previous result. 


MEASURE OF CENTRAL TENDENCY AND DISPERSION 


We need an economical way to characterize a set of measurements made 
on a system. Of course, one way would be to simply list all the numbers. 
However, this is inefficient and does not make it possible to compare 
different sets of measurements directly. On the other hand, all data can be 
described in terms of two parameters: the central tendency of the data, 
and the dispersion of the data about the point of central tendency. 


The measure of central tendency is defined as a point about which the 
data tends to group. This point is simply the arithmetic average of the data 
which is called the mean, There are two other measures of central tend- 
ency (the median and mode) but we will not need these here. You can 
pursue the subject in any introductory text on statistics if you are 
interested, 


There are two ways to compute a mean. They are defined in quite 
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different ways but both are very useful. If we are given a set of data x,, 
X9,X3,° °°, Xj, °° %, Xp, the arithmetic mean is defined as 


Xx =—2D Xj. (7) 
AN 


if you are not familiar with the sigma (2) notation in (7) it is merely a 
concise way to indicate that all the values of x; are to be summed. The 
mean is read as x bar.” 


As indicated above, there is a second definition of the mean that is very 
useful in statistical mechanics. Suppose that measurements of a system will 
yield a set of numbers x1, X2, X3,° °°, Xj, ° * °, Xp» Moreover, suppose we 
know that the system yields the measurements x; with probability P(x;) 
for i = 1,2,3, © © ©, n. Then, the mean value of the measurements of the 
system is 


x= 2D x {P(x;) Z (8) 


A simple example can illustrate the point that the mean does not com- 
pletely describe a set of data. Consider the two sets of data 


{49,50,51} and {0,25,50,75,100} . 


Both sets of data have a mean of 50. However, it is obvious that the sets 
are quite different. The second set is spread out about the mean much 
more than the first. What we need is a measure of this spreading out or 
dispersion. 


We might try to measure the dispersion by subtracting the mean from each 
element of data and then summing these differences which are called 
deviations. For the two data sets above we compute 


(x; - X) 
which gives 


{-1,0,1} and {-50,-25,0,25,50} . 


If the deviations in both sets are summed we see that 


Sum of Deviations = 2 (x; -x)=0., (9) 


Thus, the deviation cannot be used to measure dispersion. As a matter of 
fact, (9) is an alternative way to define the mean. 


A simple change is all that is needed to produce our desired measure of 
dispersion. Instead of summing the deviations, we sum the square of the 
deviations and then divide by the number of measurements. The result is 
called the variance. The variance is therefore the mean of the square of the 
deviations of a set of data. This is usually called the mean square deviation. 


Variance = s” -1y (x; - x)? (10) 


For reasons that we will not explore here, it turns out that (10) is a biased 
estimate of the variance. The unbiased estimate is 


Variance = s* = aie > (x; - x)? . 
(n- 1) 
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Since you may encounter this slightly different form, you should know 
the reason for the difference. However, we will use the first form given 
by (10). 


The variance can be computed using (10). However, x is required which is a 
separate calculation. Some straightforward algebra is all that is needed to 
convert (10) to a form that invoives only the raw data. 


n = ae ~ (Lx.)? 
Cea (11) 


n? 


Since we do not have to compute x, (11) is generally easier to use. 


Just as we have a relationship for x in terms of P(x), we should have the 
same type of relationship for s*. This is 


s* = © P(x;) (x; - x)? , (12) 
or the equivalent form which requires only raw data 

2 2 2 

S = 2x, P(x.) - (2 x;P(x;))* . (13) 


You may have wondered why in equations (10) through (13) the variance 
is written as s*. The reason is that the square root of the variance is the 
standard deviation which is usually designated by s. Recall that the 
variance is called the mean square deviation. Consequently, the standard 
deviation is called the root mean square deviation. 


MONTE CARLO METHOD 


The Monte Carlo Method is characterized by the utilization of random 
numbers, Many different types of phenomena can be investigated using the 
method, We are particularly interested in this method since the random 
Processes which are fundamental in statistical mechanics can be investi- 
gated using Monte Carlo techniques. 


The heart of the Monte Carlo Method is a sequence of numbers which 
appears to be drawn at random from a given distribution. Most computers 
have a random number generator which can be used to produce such a 
sequence. A BASIC program and printout to generate 10 random numbers 
is shown in Figure 1. 


106 REM FIGURE }! 
116 FOR I=! TO 18 
120 PRINT RWDCB) 
136 WEXT I 

$99 END 


RUN 
t 52 682 E-85 
o 386892 
o 388412 
1.64799 E-83 
6.1 7992 E-83 
e 322246 
0577867 
266971 
0981625 
«585415 


Figure 1 — Random Number Generation 
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Two interesting characteristics are seen in the printout. First, there seems 
to be no particular pattern in the numbers. Second, the numbers all fall 
between 0 and 1. The random number generator in BASIC is designed to 
conform to these two features. The numbers occur randomly in the range 
0 to 1. All numbers in the range have an equa! probability of occurring. In 
the program itself, the random numbers are generated in line 120. The 
argument of the RND function is taken to be zero. It is a dummy 
argument and does not affect the random numbers that are generated. 


Suppose we want to use the random number generator to simulate 
throwing a pair of dice. To do this we need a sequence of integers 
randomly chosen in the range 1 to 6. A program to simulate 10 tosses of a 
pair of dice is shown in Figure 2. 


1868 REM FIGURE 2 


116 =PRINT 
120 PRINT * DIE 1°,"DIE 2” 
158 PRINT 


148 DEF FNACX) =INTC6®RNDC(X)4+1) 
150 FOR I=1 TO 16 
160 PRINT FNA(G), FNACO) 


178 NEXT I 

999 END 

RUN 

DIE 1 DIE 2 
i 4 
4 i 
t 4 
4 2 
6 4 
3 3 
5 t 
2 2 
3 i 
3 3 


Figure 2 — Dice Simulation 


The desired set of random integers is generated in tine 140, With a little 
thought you can see the purpose of each part of the line, and can generate 
any other type of sequence desired. 
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THE EHRENFEST GAME 


At a somewhat primitive level, some of the ideas of random processes and 
their effect upon statistical systems can be illustrated with a game credited 
to the German physicist Paul Ehrenfest (1880-1933). In the simplest form, 
N particles are divided between two halves of a container. A move is made 
by selecting one of the N particles at random, and moving it to the other 
side of the container, A large number of moves are made in this fashion. If, 
initially, all the particles are in the same side of the container, the 
overwhelming probability is that the flow of particles will be toward the 
empty side. However, as the number of particles increases on the empty 
side, so does the probability that some particles will start moving in the 
reverse direction. This gives rise to the statistical scatter which is an 
important characteristic of macrosystems. 


A simple BASIC program to play the Ehrenfest game is given in Figure 3. 
We will use this program in some of the exercises, so it is important that 
you understand its features. The program is organized to handle 100 
particles. Each particle is represented by an element of the array B. The 
array is dimensioned in line 1020, and has all its elements set equal to zero 
in line 1030. The program segment from lines 1110 to 1150 starts the 
needed sequence of random numbers at a different point each time the 
program is run. This is done by “‘throwing away” an arbitrary quantity of 
random numbers. The program segment from lines 1210 through 1270 
enables you to set the desired number of particles on the left side of the 
container, 


1060 REM FIGURE 3 

101@ REM EXRENFEST’S GAME 

1920 DIM BI!8e8 } 

1038 MAT B=ZER 

1108 PRINT 

$118 PRINT “INPUT INTEGER BETWEEN 1 AND 188° 
1126 INPUT N 

1136 FOR I=! TO N 

$148 LET X=RND(6) 

1150 =6—NEXT I 

1208 PRINT 

1218 PRINT “TOTAL NUMBER OF PARTICLES IS 160" 
1228 PRINT “HOW MANY DO YOU DESIRE ON LEFT?” 
1238 INPUT L 

1248 FOR I=! TO L 

1250 LET Bil }=! 


1268 NEXT I 

12708) «=LLET R=108-L 
1286 PRIWT 

1290 PRINT 


13800 FOR I:i TO 56 

1318 PRINT L3 

1328 FOR J=! TO 19 

1338 LET KsINTC16@*eRND(O)+1) 
1540 IF BIK)=! THEN 1390 
§35@ LET BK)=s1 

1368 LET L=L+! 

1378 LET R=Re-} 

1388 GOTO 1422 

1399 LET BI{k1:28 

1400 LET LeL-!} 

1418 LET R=R+} 


1428 NEXT J 
1438 NEXT I 
9999 END 


Figure 3 — Program For The Ehrenfest Game 
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Of course this should be equal to or less than 100. If a particle is on the 
left, its value in the array B is 1. If on the right, it has the value 0. Thus, 
B73 = O indicates that particle number 78 is on the right side of the 
container. 


The core of the game is contained in lines 1300 through 1430. The 
program is set to print out the number of particles on the left side 50 
times, with 10 moves to be made between each printout. The random 
nature of the process is introduced in line 1330. K is randomly assigned an 
integer value in the range 71 to 100. K is then used as a subscript to locate a 
particle (an element in the array B). The particle is tested to see which side 
it is on, then moved to the other side, and the counters indicating the 
number of particles on each side are adjusted accordingly. 


Figure 4 shows two typical printouts for the game. The first started with 
all the particles on the left side. Note the steady progress towards the 
equilibrium value of 50. The second case started with 50 particles on each 
side and illustrates the scatter about equilibrium. 


RUN 


INPUT INTEGER BETWEEN { AND 188 
798 


TOTAL NUMBER OF PARTICLES IS 188 
HOW MANY DO YOU DESIRE ON LEFT? 
7168 


106 «= 98 82 74 78 68 66 62 62 68 68 58 
52 32 32 54 32 38 48 4s 46 48 58 58 
68 36 32 36 34 38 38 32 36 36 68 36 
54 52 34 36 38 46 48 32 32 30 44 48 
48 48 


RUN 


INPUT INTEGER BETWEEN 1 AND 186 
765 


TOTAL WUMBER OF PARTICLES IS 188 
HOW MANY DO YOU DESIRE ON LEFT? 
238 


38 32 52 48 32 32 34 56 34 52 58 34 
34 38 34 48 48 38 38 38 48 46 aa 44 
46 48 32 5€ 34 32 32 38 48 38 38 38 
48 a. 48 32 36 56 56 32 52 38 32 4g 
48 


Figure 4 — Printouts For The Ehrenfest Game 


EXERCISE 1 — Flipping Coins 


Suppose that 5 fair coins are simultaneously flipped. Calculate the prob- 
ability that exactly n of the coins will be heads for all possible values of n. 
Now write a BASIC program to simulate tossing the 5 coins. Arrange for 
your program to simulate 500 tosses of the coins. Each 50 tosses print out 
the frequency of occurrence of n heads divided by the total number of 
tosses, Compare the sequence of computer results to the theoretical values 
you calculated, What kind of generalizations can you make? 


EXERCISE 2 — Dishonest Coins 


Repeat Exercise 1, but assume the coins are dishonest with P(H) = 2/3 and 
P(T) = 173. 
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EXERCISE 3 — A Data Set 


Write a program to compute the mean and variance of a set of data. 
Assume that the data is contained in DATA statements, and that the last 
piece of data is followed by the “flag value’ 9999, Use the flag to detect 
when all the data has been processed. Test your program on the following 
data. 


9.7, 10.3, 9.8, 9.9, 10.1, 


EXERCISE 4 — A Frequency Distribution 


Modify the program of Exercise 3 to compute the mean and variance of 
the data given in a frequency distribution. Test your program on the 
following distribution. 





EXERCISE 5 — A Probability Distribution 


Suppose we know that a variable x occurs with the following probability: 


/ 
2 
3 
4 
5 
6 
7 
& 
9 





If we measure a great number of values of x, what will the mean and 
variance of our measurements be? Write a BASIC program to compute and 
print out the mean and variance. 


EXERCISE 6 — A Random Walk 


Write a BASIC program to simulate a 100-step random walk, Assume that 
the walk begins at x = O, and that the probability of a step forward is the 
same as a step backward. 


EXERCISE 7 — A Biased Random Walk 


A very interesting problem has been suggested by Weinstock”. Consider a 
sequence of squares numbered O through 10, A drunk is placed on square 
N. Assume that the probability of the drunk moving to a higher-numbered 
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square at each step is 0.6, and the probability of his moving to a 
lower-numbered square is 0.4, What should N be such that it is equally 
likely that the drunk will step off either end of the sequence of squares? 
Write a program to simulate many walks with different starting points to 
investigate the problem. The results are surprising! 


EXERCISE 8 — The Ehrenfest Game 


The game presented in Figure 3 assumes that the particles are moved 
between two sides with equal volumes. The result is that the equilibrium 
value of particles in each side is one-half the total. Rewrite the program to 
account for unequal volumes, Test your program to ensure that the proper 
equilibrium number of particles in each side is obtained. 


EXERCISE 9 — Relative Fluctuations 


By combining programs already written, we can investigate the fluctuation 
of the number of particles about the equilibrium point in the Ehrenfest 
game. In particular, we are interested in s/x, or the ratio of the standard 
deviation to the mean. Write a program which will play the Ehrenfest game 
and output this quantity. Arrange your program so that you can compute 
S/X for various total numbers of particles. What should happen to this 
quantity as the number of particles increases? Check your answer on the 
computer. 





CHAPTER THREE: SIMULATED 
CRYSTAL STRUCTURE 


A COMPUTER GAME 


A simple game has been described by Black! which is very useful in the 
study of statistical mechanics. Let us first examine the game in a very 
simple form. Then the connection to microphysics will be established. 


Suppose we have an ordinary checkerboard with 8 squares on a side. We 
distribute N pennies on the board in any way desired. If nj j is the number 
of pennies on the square with row number i, and column number j, then it 
follows that 


yD ms N .. (14) 


The rules of the game are as follows: 


a) Select a square at random and call it A. If square A has at least 1 
penny on it, proceed to b). Otherwise, select another square at 
random. 


b) Select a second square at random and call this square B. If B turns 
out to be the same square as A, select another square at random. 


c) Move a penny from square A to B. Repeat the entire process 
beginning with a). 


EXERCISE 10 — A Checkerboard Game 


Write a BASIC program to generate a table of integers randomly chosen in 
the range 7 to 8 inclusive. Begin at some point in the table and select two 
integers to locate square A, and the next two integers to locate square B. 
As the game is played, continue to use the number table to randomly 
locate squares on the checkerboard. Distribute 64 pennies on the checker- 
board in any way you desire. Play the game for at least 100 moves. See if 
you can make any generalizations about what is happening. Try different 
starting distributions to see if this has any effect upon your general- 
izations. 


By now you have certainly concluded that playing the game by hand will 
limit the amount of information we can gain. We should also Point out 
that there is nothing particularly important about the 8 by 8 array upon 
which we are playing the game. It could be played upon any size array. 


As advertised, we will now make the connection between the game and 
statistical mechanics. The game is a simulation of energy transfer in a 
crystal structure. In this simulation, we have a two-dimensional crystal, 
but could equally well develop a three-dimensional model if desired. Each 
square represents an atom at a point in the crystal. We assume that the 
atom can be described by a harmonic oscillator that is quantized, that is, 
the oscillator can absorb or give off energy only in fixed amounts called 
quanta. The movement of pennies on the checkerboard thus corresponds 
to transfer of quanta between harmonic oscillators in the simulated 
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crystal. The number of pennies on a square corresponds to the total energy 
of the oscillator at that point in the crystal. No pennies corresponds to an 
energy of hy /2, the lowest energy possible for a quantized oscillator. 


A computer program to play the shuffling game is given in Figure 5. The 
logic of the overall program is shown in Figure 5({a). The subroutines 
which are needed are contained in Figures 5(b) through 5(d}. Since we will 
be using this program extensively, we will discuss it in some detail to 
permit you to understand the purpose of each segment and the relation to 
the game. 





1086 REM FIGURE 5 

1018 REM A QUANTUM SHUFFLING GAME 
1926 DIM B(8,81,F(358) 

1030 MAT READ B 

1048 REM CARRY OUT TALLY 

1858 GOSUB 2660 

1196 REM RANDOMIZE RANDOM NUMBER GENERATOR 
1818 PRINT 

$126 PRINT “INPUT ANY INTEGER BETWEEN 1 AND 198" 
11:38 INPUT A 

1149 FOR I=1 TOA 

1158 LET X=RNDC@) 

1160 WEXT I 

120@ REM ENTER GAME LOOPS 

1218 FOR M=1 To 6 

1220 REM PRINT BOARD AND TALLY 
1238 GOSUB 3006 

1248 FOR N=1 TO 580 

1230 REM MAKE A MOVE 

1268 GOSUB 4900 

1278 WEXT WN 

1288 REM CARRY OUT TALLY 

1298 GOSUB 28080 

1388 NEXT ™ 

1318 REM END OF GANE LOOPS 
1328 STOP 

SOO! DATA 1,8,1,1,1,!,!,! 

8982 DATA 1,1,1,8,1,1,1,1 

SBGS DATA I,tglodatol all 

8204 DATA 1,1,1,1,!8,1,1,1 

S005 DATA 1,1,1,1,28,8,8,! 

8886 DATA Lotpdptni,d,t,! 

S007 DATA I,lalstslptetst 

8008 DATA f,igtp,t,tp,lpt,! 

9999 END 


Figure 5(a) — Program For Crystal Simulation 


The array upon which the game will be played is called B and is dimen- 
sioned in line 1020. F, the array which will be used to contain the 
frequency distribution of quanta, is also set up here. The initial distribution 
of quanta on the array B is read in from the DATA statements in line 
1030. Note that to change the initial distribution of quanta we have only 
to modify the numbers in the DATA statements in lines 8001 through 
8008. 


The TALLY subroutine branched to in line 1050 is contained in Figure 
5(b). The purpose of this is to generate a frequency distribution of the 
number of squares containing 0, 1, 2, * * * quanta. This is done by using 
the number of quanta on each square as a subscript to increment an 
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element of the array F. One is added to the subscript so that we can 
handie 0 as a subscript in the array. 


2068 REM TALLY SUBROUTINE 
2018 MAT F=ZER 

2828 FOR I=! To 8 

2038 FOR Js! To 8 

2048 =LET FIB(I, J}+! )=FIB(1, 31 1 
2850 WEXT J 

2068 NEXT I 

2678 REM END OF TALLY SUBROUTINE 
2888 RETURN 


Figure 5(b) — TALLY Subroutine 


The main loop to control the game is opened in line 1210 and closed in line Note: Lines 1100 through 1160 are needed 
1300. The strategy is to print out the array and frequency distribution of only for HP Educational BASIC and 
quanta, make 500 moves, then repeat this cycle 5 times. The OUTPUT sub- 9830 systems. 


routine is in Figure 5(c) and can be followed without difficulty. 


5888 REM OUTPUT SUBROUTINE 
5018 PRINT 
5628 PRINT 
5822 PRINT 
5824 PRINT 
SOS6 PRINT “AFTER “3(M-1)*5083"MOVES THE RESULT IS" 
S840 PRINT 
5838 FOR I=! TO 8 
50GB 4=6FOR Js! To 8 
3876 PRIWT B(I,J); 
NEXT J 


S@98 PRINT 


5130 PRINT 

5168 PRINT “QUANTA PER” ,° FREQUENCY” 
3178 PRINT “SITE 

3188 PRINT 

5196 FOR I-13 To 30 

$208 IF Fil}=8 THEN 3228 

S216 PRINT 1-3, F{1) 


NEXT I 
3238 REM END OF OUTPUT SUBROUTINE 
5248 RETURN 


Figure 5(c) — OUTPUT Subroutine 


The heart of the game is the subroutine for moving quanta which is given 
in Figure 5(d). 


4000 REM SUBROUTINE FOR QUANTA MOVES 
4010 LET IsINTC8*RNDCO)+41) 
4826 LET J=sI NTCS8*RND( 84+!) 
4050 LET K=I NTC8*RNDCO)+! ) 
4048 LET L=INTC8*RND(O)41) 
4058 IF BII,J}=8 THEN 4016 
4060 IF I#K THEN 4898 
4076 IF J#L THEN 4890 

4680 GOTO 4618 

4699 LET BI, JI=sBlI,Jl-! 
4108 LET B(K,LI=BIK,L 1 
4116 REM END OF SUBROUTINE 
4126 RETURN 


Figure 5(d) - MOVES Subroutine 
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You should study this subroutine carefully to see that it in fact follows the 
rules for the game which have been previously set down. Typical output 
from the computer game is shown in Figure 6. Only the first and last 
printouts are given. The first printout gives the initial distribution of 
quanta. The second gives the distribution after 2500 moves. 


RUN 


INPUT ANY INTEGER BETWEEN 1! AND 18 
785 


AFTER 6 MOVES THE RESULT IS 


j ! I i ! i i 1 
I I i i i ] i 1 
1 ! i i i i I i 
j i I i i t i i 
t ! i ! i ! i I 
i a 3 i i i J ! 
I 1 i i i i i 1 
i ! i ! i J i i 
QUANTA PER FREQUENCY 
SITE 
t 64 


Figure 6(a) — !nitial Printout of Crystal Simulation Program 


AFTER 2500 MOVES THE RESULT IS 


@ @ 8 8 i g 3 3 
Q 2 8 4 3 5 @ j 
2 @ 6 é 2 i g i 
@ j i 3 1 2 6 8 
G e 0 6 i 5 2 i 
g i i Q i 8 g 3 
@ i o g @ @ 3 @ 
o @ @ j i 6 } 
QUANTA PER FREQUENCY 
SITE . 
g 33 
! 16 
2 5 
 § 6 
4 i 
5 2 
6 ! 


Figure 6(b) — Final Printout of Crystal Simulation Program 





Note carefully that the printout will be different for each initialization of 
the random number generator. (See note on page 17.) Therefore, you will 
get different numbers each time the program is run. We can, however, see 
something important in Figure 6. It is clear that the final distribution of 
quanta after 2500 moves is completely different from the initial distribu- 
tion. Also, we should note that the only rule for moving quanta is blind 
chance. 


EXERCISE 11 — A Computer Game 


Enter the program in Figure 5 into your computer and run it. Examine the 
results carefully to see what is taking place. Does the initial distribution of 
quanta seem to be a reasonable one? Suppose the game were to be put 
through several more cycles. What changes do you believe would take 
place in the output? | 


INITIAL AND FINAL ENERGY DISTRIBUTION 


The results of Exercise 11 illustrate several characteristics. First, a dra- 
matic change in the distribution of energy has taken place. We started with 
1 quantum of energy at each site, and finished with something quite 
different! The last printouts indicate that quanta are moving about on the 
array but the distribution of quanta seems to be remaining about the same. 
This distribution is such that it is most likely to find O quanta at a given 
site, and least likely to find many quanta at a given site. The frequency 
distribution of quanta per site versus number of sites appears to be a 
decreasing function of some type. Later we wil! look specifically at the 
form of this function. 


A very important question remains unanswered. Does our initial distribu- 
tion of quanta determine the final distribution? If we change the initial 
distribution will we get a new final distribution? 


The computer permits us to investigate this easily. Recall that we have 
only to modify the numbers in the DATA statements in Figure 5(a) to 
change the initial distribution. Let us agree for the time being to always 
distribute 64 quanta over the 64 sites. Later we will relax this restriction. 
At this point, we are interested solely in investigating the effect of the 
initial distribution of quanta upon the final distribution. If you desire to 
suppress the printout of the array, insert the following statement in the 
program 


3035 GO TO 3090 


EXERCISE 12 — Energy Distribution 1 


Run the program in Figure 5 with DATA statements modified to put 2 
quanta on every other site. Compare the results with those from Exercise 
17. 
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EXERCISE 13 — Energy Distribution 2 


Repeat Exercise 12 except put 8 quanta on each site along one edge of the 
array. 


EXERCISE 14 — Energy Distribution 3 


Repeat Exercise 12 except put 16 quanta on each of the four corner sites 
in the array. 


EXERCISE 15 — Arbitrary Distribution 


Repeat Exercise 12 with a distribution of quanta of your choice. If you 
put a large number of quanta on one oscillator you may have to change 
the DIM statement in line 1020 of Figure 5{(a). 


DISTRIBUTION OF ENERGY 


The computer results thus far point up a very important characteristic of 
our simulated crystal structure. No matter how the quanta are distributed 
initially, we a/ways obtain the same distribution of quanta per site after a 
very large number of moves. There is some question about just what this 
distribution is since the computer results do not settle out on some fixed 
set of values but instead tend to oscillate statistically. Our task in this 
section is to learn more about this unknown distribution. 


Suppose that the computer program to simulate random quanta inter- 


change has been running for some time. We extract one frequency distri- 
bution of quanta per site and suppose it turns out to be that shown in 


Figure 7. 





Figure 7 —- Typical Distribution of Quanta In Simutated Crystal 


Our experience with the computer simulation tells us that the outputs 
made both prior to and after the printout in Figure 7 will differ from that 
above. What we want to find is that distribution about which the indi- 
vidual distributions are scattered. If we were to take many individual 
distributions and average them, reason tells us that the result should be 
close to the theoretical! distribution. 


In any event, the distribution in Figure 7 points out rather clearly the type 
of mathematical function needed. A decreasing exponential would be an 





ideal function. Let us assume that the distribution is given by 


n= ne ibe (15) 


where ng and 8 are unknown constants, j is an integer giving the number of 
quanta per site, and € is the energy per quantum. In our exercises thus far, 
we have tacitly assumed that € = 1, although there is no need for this 
assumption. As we have set up the simulation, we generally will know N, 
M, and e. N is the total number of sites, M is the number of quanta to be 
placed on the sites, and € is the energy of each quanta. The total energy 
on the simulated crystal is given by 


E,=Me . (16) 


With this information, we want to find ng and 6 which will enable us to 
use (15) to compute the theoretical distribution. 


To begin, we note that 
N=2 ni, (17) 
or 


N=1No + noe PE + noe 2PE pote e 
=No (1 +e PE + 9 2he 4. ce) (18) 
The binomial theorem states that 
(1 +e PE yo Bey... y= (1 ~ @7Pe)-1 (19) 
provided that 
(eBE\2 < 4 


which is true provided that both 8 and ¢€ are greater than zero. We will 
assume this is the case. 


Thus 


N=ng(1- e| Be)-1 (20) 


The total energy of the system is given by 


E, = (no) (0) + (my) (e) + (ng) (2e) ++ + - (21) 


a = ne ibe ' 
which leads to 
E, = noe Ke + noe” 2€2¢ acess 


= noee PE (1 + 207 BE + Ze" 2BE 4. . wy (22) 
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Again, the binomial theorem can be used to simplify the quantity inside 
parentheses on the right. The result is 


E, = noee PE (1 - e Pe)-2 (23) 
The average energy per oscillator is E,/N, or 


€ 





E = . (24) 
QE _ 4 
Solving for B we obtain 
1 = 
B =o int +eé/E) . (25) 
If we solve (20) for ng we have 
Nog =N(1- e~ Pe) (26) 


Equations (25) and (26) enable us to derive a theoretical distribution given 
the energy per quantum e, the total number of quanta M, and the total 
number of oscillators N. For the exercises we have done thus far, € = 1, N 
= 64, and M = 64. Consequently, E = 1. Equations (25) and (26) yield 8 = 
In(2), and ng = 32. The distribution is therefore 


a9 Qj |nl2) a5 (1)! 
nj 32e 32(4) 


Figure 8 shows the comparison between our derived distribution and the 
experimental one from Figure 7. 
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Figure 8 — Comparison of Typical Computer Results 
and Derived Distribution 


We have obtained a distribution of quanta per site which agrees well with 
Figure 7, and reasonably well with other computer results. Note carefully, 





however, that everything in the devefopment flows from the assumption 
that 


sr noe ie (27) 


The fact that (27) is the correct distribution function has not been proven. 
However, it turns out that it is the correct relation. The derivation and 
proof are given in most statistical physics texts. 


Before going on, we should pause to consider the restrictions on the 
determination of 8 and ng by Equations (25) and (26). First, we used a 
binomial series to simplify expressions in two cases. This implies an 
infinite number of terms in the series. Actually, we have only from 6 to 10 
terms in the actual computer results. Therefore, the fewer the number of 
terms, the greater the disagreement between theoretical and computer 
results. Also, in the development of (20) we require that 6 and € be greater 
than zero. This has more significance than would appear at first glance and 
will be discussed later. For the present, we assume only that 6 and € are 
positive constants. 


The last item to be discussed in this section is the connection between the 
distribution function given by (27) and the probability of observing a site 
with j quanta. From (1), 

n 


n: 
i a ~ j Be 
P(j) NN e 


However, je = Ej the energy on an oscillator with j quanta. Thus 


P(j) = CePEj (28) 
where C is a constant (explained later in this discussion). The probability 
relation given by (28) is much more general than we might believe. In our 
treatment, J is an integer. More generally, if r is a specified state and E, is 


the energy associated with that state, then 


P(r) = Ce~BEr - (29) 


This expression has fundamental importance in statistical mechanics. The 


factor e BE; is known as the Bo/tzmann factor. The entire distribution 
given by (29) is called the canonical distribution. The constant C can be 
determined (in principle at least) by the condition that 


> P(r)=1. 
: 


Thus 


c=(z eBEr)-1 (30) 


EXERCISE 16 — An Average of Distributions 


Madifty the crystal simulation program to average 20 individual distribu- 
tions of quanta per site. Start the program with DATA statements which 
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reflect a typical state had the program been running for some time. 
Compare the average of the frequency distributions to the theoretical 
distribution. Can you make any generalizations? 


EXERCISE 17 — A Larger Simulation 


Modify the crystal simulation program to handle an N by WN array with M 
quanta distributed on the array. Suppress the printout of the array itself. 
Output only the frequency distribution of quanta per site. Test your 
program on a 20 by 20 array containing 400 quanta. Compute the 
theoretical distribution and compare to your computer results. 


EXERCISE 18 — Effect of Total Energy 


Use the program from Exercise 17 with a 10 by 10 array. Run the program 
with 100, 200, 400, and 1000 quanta. Examine your results carefully, and 
in each case compare to the corresponding theoretical distribution. 


NEW RULES — NEW DISTRIBUTIONS? 


We have already investigated thoroughly the idea that the initial distribu- 
tion of quanta has no effect upon the final distribution obtained after 
randomly moving quanta many times. In all cases, the exponential distri- 
bution obtained is quite independent of the initial distribution. 


An equally important question has not been answered. Is it possible that 
the rules for moving quanta determine the final distribution? Certainly, if 
we think carefully about the simulation there are several aspects which 
might seem unreasonable and could suggest new rules. First, we have 
assumed that any oscillator can transfer a quantum of energy to any other 
oscillator in the array. How does this take place? We have avoided any 
discussion of the mechanism by which energy is transferred. However, 
whatever the mechanism is, it would seem reasonable that a quantum 
should be transferred to an adjacent site with a higher probability than a 
more distant site. Also, we have been transferring energy in discrete 
amounts. What would happen if we transferred fractional parts of the 
energy on a given oscillator? 


Remember that the “‘rules’’ of the game are contained in the subroutine in 
Figure 5(d). If we change the rules, this is the primary point in the 
program that must be modified. If non-integer amounts of energy are 
transferred, then the amount of energy on each oscillator will, in general, 
be a non-integer value. We can still use the TALLY subroutine in Figure 
5(b). However, BASIC will round non-integer values to the nearest integer 
when used as a subscript. 


EXERCISE 19 — New Rules 1 


Modify the rules of the crystal simulation as follows: If the donor 
oscillator is not on the edges of the array, select the oscillator to receive 
energy at random from the four closest neighbors. If the donor oscillator is 
on the edges of the array, but not one of the four corners, select the 





oscillator to receive energy at random from the three closest neighbors. If 
the donor oscillator is a corner oscillator, select the oscillator to receive 
energy at random from the two closest neighbors. With these new rules, 
run the simulation for various-sized arrays and total number of quanta. In 
each case, compare the computer results with the Boltzmann distribution 
given by (27). Did the rule change affect the final distribution? 


EXERCISE 20 — New Rules 2 


Modify the rules of the crystal simulation as follows: When the donor 
oscillator is selected, instead of transferring a quantum of energy, transfer 
a fractional part of the energy on the oscillator. Determine randomly what 
this fraction will be. Run the program for various-sized arrays and initial 
amounts of energy. In each case, compare the computer results with the 
Boltzman distribution given by (27). Did the rule change affect the final 
distribution? 


EXERCISE 21 — Goodness Of Fit 


One of the difficulties with Exercises 19 and 20 is how to detect if there is 
a significant difference between two distributions. Consult an introductory 
statistics text and learn how to perform a Chi-Squared test to accomplish 
this. 
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CHAPTER FOUR: TEMPERATURE 
AND HEAT FLOW 


In the crystal simulation which has been the central part of our investiga- 
tion so far, we have not mentioned two important terms. The first is 
temperature, which is commonly used to describe macrosystems. The 
second is heat, which has central importance in thermodynamics. We can 
look at both these concepts using our crystal simulation. 


DEFINITION OF TEMPERATURE 


Recall that the Boltzmann distribution Is 
ni = nge ie : (31) 


and that the average energy per oscillator is 


€ 


E= 
ePE _ 





(32) 


Thus far we have only defined B as a positive constant (for a given system) 
and have not associated it with any macroparameter. Suppose we compare 
two systems that are identical except one has a higher average energy per 
oscillator than the second: 


— E€ 
ie eae, (33 
1 Bie _ 4 ) 
= € 
E, =———_-.. (34) 
ere _ 4 


lf E, > E, then it follows from (33) and (34) that B, > B,. In other 
words, as E increases, 8 decreases. Increasing the average energy per 
oscillator in the system corresponds to increasing the temperature of the 
system. This in turn must correspond to a decrease in 8. It turns out that 


1 
Bar : (35) 


The constant of proportionality in (35) is I/k where k is the Bo/tzmann 
constant. Therefore 


geek. (36) 


which connects temperature to the Boltzmann factor and the canonical 
distribution in (29). 8 must have dimensions of energy ‘ and this is indeed 
so by (36). Also recall we have specified that 8 must always be positive. 
Since k is a positive physical constant, then T must always be positive. 
Thus T as determined by (36) is the abso/ute temperature. 
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Since we can compute B from information about the crystal simulation we 
can, in turn, compute the temperature of the system, 


T ee (37) 
k In(1 + €/E) 

However, we must be careful about using (37) to compute the tempera- 
ture. The relationship is valid on/y when the system is described by the 
Boltzmann distribution. If, for example, all the quanta are on the same 
oscillator, E is the same as when the quanta are distributed according to 
the Boltzmann distribution. However, the two situations are quite dif- 
ferent. When the system is in equilibrium it is described by the Boltzmann 
distribution and the temperature can be determined from (37). 


To see the effect of temperature, suppose we have a system composed of 
elements that are known to be in one of two energy levels. These energies 
are O and € respectively. By (29) 


P(0) = Ce“ P() 


and 


P(e) = CeBE 
The constant C is determined by (30) 

C= (e°Bl0) 4 Be)” = (1 + @ Be) 
Thus 


p(a)= (1 +e7Fe) (38) 


and 


pte) =e7F (1 + 0B) (39) 


We can compute the average energy of each element in the system by the 
method described by (8). 


E = (O)P(O) + (e)Ple) , 


or 


ee” Be 


E = (40) 


1 + @ BE | 


It is a simple mathematical exercise to show that if the temperature 
approaches absolute zero, then kT ~ Q, and E — O. If the temperature 
increases without limit, kKT—> °°, and E > e/2. If there are N elements in the 
system, then the total energy of the system is simply N E. 


With this simple case there is no need to involve the computer. However, 
with more involved (and more interesting) cases where the energy level 
structure is more complicated, the computations are much more involved. 
Applications of the computer to this type of problem have been described 
by Weinstock.° 





EXERCISE 22 — Unequally-Spaced Three Level System 


A system is composed of units each of which has the permissible energy 
levels below, In each case, write a program to compute the average energy 
per unit versus kT/e. Consider at least the range 0.05 < kT/e < 3.0. 
Plot the results. 


a) 0, €, 3e 
b) 0, 2€, 3e 
c) 0. €, 5€ 


d) O, 4e, 5€ 


EXERCISE 23 — Equally-Spaced Level System 


The harmonic oscillators of our crystal simulation have equally spaced 
energy levels: O, €, 2€, 3e,°* +. Find the average energy per oscillator 
versus kT/e for the cases below. Consider at least the range 0.05 <kT/e< 
3.0. 


a) 0,€ 2€ 
b) 0, €,2€,° °°, 10€ 


c) 0.¢€,2€,°°*, 30¢€ 


EXERCISE 24 — Inverse Square-Spaced Level System 


Some systems have energy levels of the form 


C 
rae 


n=1,2,3°°"* 


in our analysis we want to have the lowest energy level be 0. The 
expression above modified to do this, and, with an assumed value of C, is 


— ] — ¢* o® @® 
E,, = 5€ ('-5) n=1,2,3, 


Write a program to find the average energy per unit versus kT/e for such a 
system. Use n = 1, 2,3, °° +, 10anda reasonable range of kT/e. Plot the 
results. 


EXERCISE 25 — Randomly-Spaced Level System 


One of the interesting things about investigations with the computer is 
that we can simulate any type of universe we desire. [t makes no dif- 
ference whether there is a real counterpoint to the simulation or not. 
Suppose that natural law was such that a system had units which had 
randomly-spaced energy levels. How would the average energy per unit 
versus kT/e behave for such a system? Assume that five energy levels are 
selected at random from the range O to 5e. 
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HEAT CAPACITY 


With the foundation established in the previous section, we can very 
quickly bring out the concept of heat capacity. Suppose that a system has 
average energy per unit of E. In the crystal simulation, the unit was an 
oscillator. In other systems, the unit might be something different. The 
essential idea is that the system is composed of units, each of which can 
hold energy, and has an average energy. The total energy of the system is 
the sum of the average values of all the units in the system. 


Exercises 22 through 25 were concerned with finding E as a function of 
kT/e. In other words, we found how the average energy per unit depended 
upon the temperature of the system. Recall that for the very simple 
system in which the units could be only in the levels 0 or e€, the average 
energy per unit was 


= -Be 
E = ee ae: , (41 ) 
1+e- Be 
where 
1 
Bray 


The heat capacity (at constant volume) per unit is defined as the rate of 
change of E with respect to T, 


C, = (3E/dT), . (42) 


Since (41) depends explicitly on 6, we can use the chain rule of differenti- 
ation to give 


_ (28) (26. 
a (ie), (3) _ 


Carrying out this differentiation and simplifying we have the following 
expression for the heat capacity per unit for the system described by (41), 


e~Be 


——— . (44) 
(1 +e Pey? 


C= k (Be)? 


It will simplify our computer exercises if we compute C,/k instead of C,. 
Thus, 


C - Be 
_V_ es 


For high temperatures, Be << 1 and C,)/k falls off as the inverse square of 
the temperature. For low temperatures, Be >> 1 and C,/k rises expo- 
nentially. 


The only reason we can write (45) is that the functional dependence of E 
upon B (and therefore T) is known. Generally this information will not be 
available to us. However, we can compute and tabulate E versus kt/e. 
Using the computer we can then numerically compute the derivative of E 





with respect to T: 








o -GE__0E_ dlkT/e) _ k _ 8 
VOT dO(kT/e) dT € 0(kT/e) 
or, 
e. _ dE 
k “vy > BikT/e) (46) 


Now we need to compute numerically the right side of (46). Suppose that 
kT/e takes on the values .05, .10,° °° , 2.95, 3.00. Ifi=1,2,3,° °° , 60, 
then the values of kT/e are simply .051. We can compute. E; for each of the 
values of i. The first derivative of E with respect to kT/e (using a central 
difference expression) is2: 

oe; Ein Gi-y 
d(kT/e)  2A(kT/e) 


But since A(kt/e) = .05, we have 


Cy\ | Fita > Fi-t ne 
k /; 1 | 


We have dropped the € on the left since E is determined in terms of € and 
they cancel. Using (47) we can compute C\/k at each of the points in the 
array of kT/e values to get a profile. Remember that we must multiply the 
computed values by k to get the heat capacity per unit, and multiply by 
Nk to get the heat capacity of the system, where N is the number of units 
in the system. 


The primary purpose here is not to compute numerical values of heat 
capacity. Rather, it is to compute the draw profiles of how the heat 
capacity varies with temperature and thereby gain insight into the behavior 
of the systems. 

EXERCISE 26 — Heat Capacity, Unequally-Spaced Three Level System 
Write a program to compute the heat capacity per unit for each of the 
systems in Exercise 22. Plot your results. 

EXERCISE 27 — Heat Capacity, Equally-Spaced Level System 

Write a program to compute the heat capacity per unit for each of the 
systems in Exercise 23. Plot your results. 


EXERCISE 28 — Heat Capacity, Inverse Square-Spaced Level System 


Write a program to compute the heat capacity per unit for the system in 
Exercise 24. Plot your results. 


EXERCISE 29 — Heat Capacity, Randomly Spaced Level System 


Write a program to compute the heat capacity per unit for the system in 
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Exercise 25. Plot your results. 


HEAT TRANSFER 


The crystal model furnishes an ideal way to study heat transfer. First, we 
must be very specific in our definition of heat. Heat is energy transferred 
in the absence of any external work. In our crystal model, we could 
arbitrarily divide the system (the crystal) into two smaller subsystems. We 
could take these to be the left and right halves of the crystal respectively. 
Now, as quanta are transferred between the oscillators, energy is clearly 
flowing from one subsystem to the other. It is also clear that no external 
work is involved. As indicated above, the energy thus transferred ts defined 
as heat. 


Several important questions can be investigated using the simulated crystal 
model. Suppose two systems at the same temperature are brought together 
and allowed to come to equilibrium. Will the distribution of quanta differ 
from distributions in the two separated systems? What will the outcome be 
if the two different systems are initially at different temperatures? Finally, 
suppose a small system at equilibrium and at one temperature is brought 
into contact with a large system at equilibrium at another temperature. 
What will take place? 


By this time you should be convinced that if we place M quanta on N 
oscillators and move them about randomly, after a sufficiently-long time 
the quanta will be distributed according to the Boltzmann distribution and 
we can compute the temperature using (37). In the exercises to follow, the 
computer is not needed. You can reach a solution to each of the exercises 
using your knowledge of how blind chance affects systems, the fact that E 
is the average energy per oscillator, and equation (37). 


EXERCISE 30 — Two Systems, Same Temperature 


Two systems at the same temperature are isolated from one another for a 
long time. They are then brought together and allowed to reach equilib- 
rium, What will the temperature of the combined system be? Will any heat 
flow take place? 


EXERCISE 31 — Two Systems, Different Temperature 


System A is at equilibrium and has temperature T A System B is also at 
equilibrium and at temperature T B Suppose that Th > T If the two 
systems are brought into contact and permitted to exchange energy, what 
will the final temperature be? Will heat flow? If so, in which direction? 


EXERCISE 32 — Two Systems: One Large, One Small 


A small system at equilibrium is isolated from a large system also at 
equilibrium. Suppose that the temperature of the small system is not the 
same as that of the large one. If the two systems are brought together, 
describe what will take place. 





CHAPTER FIVE: ENTROPY AND 
EQUILIBRIUM 


ACCESSIBLE STATES AND ENTROPY 


The introduction pointed out that the bridge between macrostates and 
microstates is furnished by the relation 


S = kIn(QQ) , (48) 


where S is the entropy of a system, k is the Boltzmann constant and {2 is 
the number of microstates accessible to the system. The new concept here 
is that of accessible microstates. An example will illustrate the essential 
ideas: Suppose we have a system composed of three oscillators which share 
a total of 5 quanta, where each quantum of energy has magnitude e. 
Figure 9 enumerates all possible combinations which are allowable. 


Configuration 





Figure 9 — Configurations and Microstates for Three Oscillators 
Sharing Five Quanta of Energy 
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There are five distinguishable configurations for the system, and a total of 
21 allowed microstates. The configurations are determined by the 
frequency distribution of quanta per site. These distributions are shown in 
Figure 10. 













0 Z 
5 1 
0 1 
1 1 
4 1 
0 1 
2 1 
3 1 
1 2 
3 1 
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Figure 10 — Distribution of Quanta for Possible Configurations 


The frequency distribution of quanta per site identifies each of the 
possible configurations. This is just the familiar frequency distributions 
which were output from the crystal simulation. We can use this 
information to calculate the number of microstates corresponding to each 
configuration: 


N! 


8 ED (fal) (gl 


(49) 


in (49), N is the number of oscillators and f,, fz, fs, °° ° are the 
frequencies in the distribution of quanta per site. From Figure 10 we have 


= opm 73 
% = Gays 3: 
% = ao 73 


These computations agree exactly with the results enumerated in Figure 9. 
Using (49) we can compute Q for any configuration. However, rather than 
computing Q, we need In(Q2) since this is required to compute the entropy 


S. Taking the natura! logarithm of both sides of (49) we have 


In(&2) = In(N!} - In(f,!) - In(f5!) - In(f4!) te (50) 


It is not possible to evaluate the factorials in (50) on the computer and 
then take the required logarithms, because of the size of factorials. If, as in 
our first crystal simulation, N = 64, then N! = 1.2689 x 10®?. This very 
large number is outside the range of many computers. We can solve the 
problem by using Stirling’s approximation for n!. 


In(n!) = ( + p)Inin) ~ n+ in V2 ; 


2 


which if the numerical value for In ./ 277 is inserted becomes 
In(n!) = (r + 5)in(n -n+0.918939 . (51) 


The approximation gets better as n increases. For n= 5 the error is about 
16 parts in ten thousand. For n= 50, the error is about 2 parts in a 
million. The error for small n does not play an important part in the 
computer investigation of entropy. Consequently we will always employ 
Stirling’s approximation given by (51) to compute the logarithms of the 
factorials required. 


It will be simpler to solve for S/k rather than S. Thus 


x|” 


= In(N!) - In(f,!) - In(f5!) a. eae (52) 


where each of the terms on the right is computed using (51). 


EXERCISE 33 — Entropy Calculation 1 


Modify the program for the crystal simulation to produce a printout of 
S/k for each of the frequency distributions as they are computed. Set up a 
10 by 10 array which contains 


a} 200 quanta 


b) 100 quanta 
Examine the results carefully and try to generalize what ts taking place. 
Experiment with various numbers of moves between printouts of S/k. 
EXERCISE 34 — Entropy Calculation 2 


Repeat Exercise 33 except use a 5 by 5 array which contains 
a) 50 quanta 


b) 25 quanta 
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EQUILIBRIUM 


We are now in a position to firm up previous discussions about equilib- 
rium. First, we must return to the information in Figure 9 to make a 
fundamental point concerning statistical mechanics. We had 3 oscillators 
sharing 5 quanta of energy which gave rise to 21 allowable microstates. 
The key postulate of statistical mechanics is that all microstates are 
equally likely. Thus, configuration 1, which contains 3 of the possible 21 
microstates, will be observed with a probability of 3/21. Likewise, the 
probability of observing the system in configuration 3 is 6/21, and so on. 
It follows that the most probable configuration is the one corresponding 
to the largest number of microstates. 


Suppose we look at a typical example to make this point clear. Two 


frequency distributions of quanta per site are given below. In each case, we 
will assume N = 64, and 64 quanta. 


Quanta Per Quanta Per 


16 
8 
5 
3 
1 











From (49), we have 


«gal 
2, = (63) 4h ah 64, and 

vy ee Ot 32 
9 = Bonen 8) en) enayn (22 x 19 


Therefore, the probability of configuration 2 !s about 10°} greater than 
the probability of configuration 1. The system will be found with over- 
whelming probability in that configuration which corresponds to the 
greatest number of microstates. It is possible but not probable to find a 
system in a configuration with a small number of microstates. Just how 
small this probability is can be seen in the computations for Q, and Q, 
above. 


The equilibrium configuration of a system is that one which corresponds 
to the greatest number of microstates. The Boltzmann distribution of 
quanta per site corresponds to a greater number of microstates than any 
other distribution. if a system is permitted to respond to blind chance in 
the interchange of quanta, it will drift towards the equilibrium configura- 
tion which implies the largest number of microstates. 


Since S/k is the natura! logarithm of the number of accessible states, we 
can immediately conclude that as a system goes from non-equilibrium 
toward equilibrium, the entropy will increase. Moreover, the equilibrium 
configuration corresponds to a maximum in entropy for the system. 





EXERCISE 35 — Entropy and the Boltzmann Distribution 


Figure 8 contains the Boltzmann distribution for 64 quanta on 64 oscil- 
lators. The first 6 terms are: 32, 16,8, 4,2, and 1. These terms add up to 
57 quanta on 63 oscillators because of truncation error. Write a program 
to compute S/k. Use this program to investigate the effect of changing the 
distribution (keeping 57 quanta and 63 oscillators). Compare to the 
Boltzmann distribution. 
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CHAPTER SIX: FINAL THOUGHTS 


There have been a number of assumptions made throughout this unit. As 
long as the subject is restricted to classical statistical mechanics, it is not 
necessary to raise the issue. However, as soon as we consider quantum 
statistics, it is important to understand just what these assumptions are. 


We have assumed that the oscillators of our simulated crystal were identi- 
cal in all respects save that of location. We could identify each oscillator 
by position and determine the number of quanta of energy held by that 
oscillator. In other words, the system was composed of a number of 
distinguishable units. Such a system is described by Boltzmann statistics as 
we have seen. 


If, however, we consider a system of particles which are not distinguish- 
able, the situation is significantly different. An example might be a group 
of identical atoms moving inside a rigid box. The atoms move freely in a 
common potential well which is not divided up into local valleys (which 
would correspond to our oscillators in the crystal simulation}. Such 
particles cannot be assigned to a location in the box. In principle, we can 
observe the energy of every particle but there is no way to determine 
which of the particles had which of the observed energy values. Thus the 
particles are indistinguishable. 


If the set of indistinguishable particles obeys the exclusion principle (no 
two particles can be in the same dynamical state) then the system is 
described by Ferm/i-Dirac statistics. If the set of indistinguishabie particles 
is not restricted by the exclusion principle, then the system is described by 
Bose-Einstein statistics. Quantum statistics is concerned with these two 
new types of statistics. 
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ANSWERS TO SELECTED 
EXERCISES 


CHAPTER TWO 
EXERCISE 1 — Flipping Coins 


The theoretical results are P(O) = P(5) = 0.03125, P(1) = P(4) = 0.15625, 
and P(2) = P(3) = 0.31250. The program and results after 451 throws 
follow. 


LIST 

108 REM EXERCISE ! 

118 DIM F(6) 

128 MAT F=ZER 

138 PRINT 

148 LET c= 

150 FOR I2i TO 568 

168 IF C<50 THEN 248 

170 PRINT “AFTER*sI3“TOSSES THE RESULTS ARE” 
188 PRINT “HEADS”, PROBABILITY” 
190 PRINT 

2@6 FOR Jz! TO 6 

218 PRINT Jl, F(J)/I 

226 WEXT J 

225 PRINT 

226 PRINT 

258 LET C2 

248 LET Ss! 

258 FOR Ki TO 5 

268 LET SsStINTCRND(8)4+.5) 
276 WEXT XK 

289 LET F(S}2F(S}! 

296 LET C=C+l 

308 NEXT J 

999 END 


RUN 


AFTER 451 TOSSES THE RESULTS ARE 
HEADS PROBABILITY 


51842 E-82 
0195122 
277162 
0538377 
213969 

2 45982 E-82 


Wb WN=— & 


EXERCISE 3 — A Data Set 


100 REM EXERCISE 3 

118 LET NeSi:=S2=¢ 

120 READ X 

138 IF X=9999 THEN 182 

140 LET N=N+1 

158 LET Si=Si+x 

168 LET S2:=f4+xt2 

176 GOTO 1286 

188 PRINT 

198 PRINT “MEAN = “sSi/N 

262 PRINT “VARIANCE = "3 (Ne S2-S!1 #2)/Nt2 
eos at 9 of, 18 63,968,929 ,18.1,18,108.2,9999 
99 


RUN 


MEAN = 9.91429 
VARIANCE = .1383543 
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EXERCISE 5 — A Probability Distribution 


103 REM EXERCISE 5 

110 LET Si:<S23 

120 READ RW 

136 FOR Ist TO N 

148 READ X,P 

138 LET SisSi+XxX*P 

168 LET S2=S2?Xt2eP 
178 WEXT I 

188 LET $8} 

192 LET V=S2e-S! % 

200 PRINT ° MEAN = "3M 
210 PRINT * VARIANCE = °3V 
886 DATA 9 

801 DATA 1,4.80000 E-92 
8822 DATA 2,8.80000 E-82 
803 DATA 3,.12 

884 DATA 4,.!16 

605 DATA 5,2 

8@6 DATA 6,.16 

887 DATA 7,12 

808 DATA 8,8.88000 E-82 
8@9 DATA 9,4.808800 E-92 
999 END 


RUN 
MEAN = 5§ 
VARIANCE = 4 


EXERCISE 7 — A Biased Random Walk 


According to the results of three different starting points, the answer iS 
between step 1 and 2. 


100 REM EXERCISE 7 

110 PRINT "INPUT STARTING POINT"; 
128 INPUT X!1 

130 LET L=R=8 

140 FOR f-1 TO 180 

145 LET X=xl 

158 IF RND(@)<.6 THEN 198 
1780 LET X=Xe! 

186 G0TO 206 

199 LET X2X+} 

2600 IF x<@ THEN 236 

216 IF X10 THEN 259 

226 GOTO 138 

230 LET L=L+! 

248 GOTO 26 

250 LET RsR+! 

260 NEXT 2 

278 PRINT °PCL) = “31/188 
286 PRINT "PCR) = “sR/188 
298 PRINT 

366 GOTO 118 

999 EMD 


RUN 
INPUT STARTING POINT?76 


PCL = 62 
PCR) = 38 
INPUT STARTING POINT?! 
PCD = 4! 
PCR) = 259 
INPUT STARTING POINT?72 
PCL) = eZ] 
PCR) = 279 





EXERCISE 9 — Relative Fluctuations 


108 REM EXERCISE 9 

110 DIM Blise) 

120 MAT B=ZER 

130 PRIWT 

140 PRINT “INPUT INTEGER BETWEEN 1 AND !98° 
156 INPUT Z 

1628 FOR Is] TO Z 

170 LET X=RND(B) 

180 NEXT I 

190 PRINT “INPUT TOTAL NUMBER OF PARTICLES” 
200 PRINT “MUST BE EVEN AND NO MORE THAN 160" 
218 INPUT N 

220 LET LeR=N/2 

238 FOR Isi TO L 

246 LET Bt! )=} 

250 NEXT I 

268 PRINT 

276 LET Sj=S2: 

28@ FOR I=1 TO 168 

298 LET S1=Si14L 

386 LET S2:=S24+Lt2 

310 LET K=I NTC NeRND(O)+41) 

328 IF B(X)=i THEN 388 

338 LET BIK1is! 

346 LET L=L+t 

568 LET R=Re! 

378 @QOTO 418 

380 LET B(K )}=@ 

398 LET L=Le-l 


420 LET M=S1/100 

438 LET S=SQR( (100 2 -S1*S1)/19668) 

a pore “RATIO OF STND DEV TO MEAN = “3S/M 
99 D 


RUN 


INPUT INTEGER BETWEEN 1 AND 168 
267 

INPUT TOTAL NUMBER OF PARTICLES 
MUST BE EVEN AND NO MORE THAN 108 
748 


RATIO OF STND DEV TO MEAN = ,178249 


READY 
RUN 


aoe INTEGER BETWEEN ! AND 168 
? 


INPUT TOTAL NUMBER OF PARTICLES 
MUST BE EVEN AND NO MORE THAN 108 
268 


RATIO OF STND DEV TO MEAN = ,.114215 


READY 
RUN 


INPUT INTEGER BETWEEN ! AND 180 
781 

INPUT TOTAL NUMBER OF PARTICLES 
MUST BE EVEN AND NO MORE THAN 182 
788 
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EXERCISE 9 — (Continued) 
RATIO OF STND DEV TO MEAN = 6.135737E-02 


READY 
RUN 


INPUT INTEGER BETWEEN 1 AND 108 
243 

INPUT TOTAL NUMBER OF PARTICLES 
MUST BE EVEN AND NO MORE THAN 168 
2108 


RATIO OF STND DEV TO MEAN = 5,819357E-@2 
CHAPTER THREE 


EXERCISE 11 — A Computer Game 


The first and last printouts are contained in Figures 6(a) and 6(b) respec- 
tively. After 500 moves, the frequency distribution is oscillating about the 
equilibrium value. Results in this exercise and all others that utilize the 
random number generator will be different depending upon how the 
random number generator is sequenced in your computer. The trends 
should remain the same, however. 


EXERCISE 13 — Energy Distribution 2 


Modify the DATA statements in the program in Figure 5 as shown below. 
The distribution goes over to the same type as for Exercise 11. 


8801 DATA 


B_8, 

8002 DATA @,0,0,0,0,80,0,0 
8003 DATA @,8,0,0,0,0,0,0 
8004 DATA 8,0,0,0,0,9,0,2 
8005 DATA 0,0,8,0,0,0,0,0 
8806 DATA @,0,0,9,0,0,0,0 
8807 DATA 9,8,0,0,0,0,0,0 
8808 DATA 2,2,0,0,0,0,0,0 
9999 END 


EXERCISE 15 — Arbitrary Distribution 


The guanta can be distributed in any fashion desired on the oscillators. 
After a sufficiently large number of moves, all initial distributions will 
drift over to the type of distribution seen in Exercise 11. 


EXERCISE 17 —A Larger Simulation 


The complete program is given below. Depending upon the size of the 
program space available to you, the DIM statement may have to be 
modified. Also, if you have sufficient program space you can continue to 
specify the initial configuration with DATA statements. With limited 
space, other provisions have to be made. Following the program, the first 
few printouts are shown. The printouts are made every 50 moves to watch . 
the process more closely. 


1000 REM EXERCISE 17 

1918 REM LARGER SIMULATION 

1026 REM MAX SIZE IS 20 BY 26 

1038 DIM B(20,28),F{(5d) 

1940 PRINT “INPUT SIDE OF SQUARE ARRAY"; 
165@ INPUT Ni 





EXERCISE 17 — (Continued) 


1660 
1078 
1682 
1188 
1118 
1128 
1138 
1148 
1156 
1168 
1290 
1218 
1220 
1238 
12408 
12528 
1268 
1278 
1280 
1296 
1388 
13198 
15268 


2080 
2018 
2828 
2238 
20648 
28 58 
28608 
2078 
2888 
5888 


38268 
5658 
3848 


3860 
5878 
5688 
35890 
5188 
S116 
5128 
31 56 
51 40 
46006 
4018 
4628 
4230 
48 48 
40 56 
4268 
49708 
4086 
4096 
4108 
4110 
4120 
9999 


MAT BoCON[NI,N1) 

REM CARRY OUT TALLY 

GOSUB 20600 

sa ane RANDOM NUMBER GENERATOR 
N 

PRINT “INPUT ANY INTEGER BETWEEN 1 AND 160"; 

INPUT A 

FOR I=! TO A 

LET X=RNDCO) 

NEXT I 

REM ENTER GAME LOOPS 

FOR Miz! TO 122 

REM PRINT TALLY 

GOSUB 3008 

FOR N=i TO 58 

REM MAKE A MOVE 

GOSUB 4896 

NEXT N 

REM CARRY OUT TALLY 

GOSUB 2002 

NEXT Mi! 

REM END OF GAME LOOP 

STOP 


REM TELLY SUBROUTINE 

MAT F=Z ER 

FOR Ii TO Ni 

FOR J=i TO NI 

LET F(B({I,J#! EFIB(I, JH! #1 
NEXT J 


NEXT I 

REM END OF TALLY SUBROUTINE 
RETURN 

REM OUTPUT SUBROUTINE 

PRINT 


PRINT 


PRINT 

ee “AFTER” 3 ( Ml-1)*503;”" MOVES THE RESULT IS” 
PRIN 

PRINT “QUANTA PER” ,” FREQUENCY” 
PRINT “SITE” 

PRINT 

FOR I=1 TO 38 

IF F(I1}=6 THEN 3126 

PRINT I-1, F(T) 

NEXT I 

REM END OF OUTPUT SUBROUTINE 
RETURN 

REM SUBROUTINE FOR MOVES 

LET ISINTCNI*RNDCO) +1) 

LET J=I NTC NI*RNDCOD+8) 

LET K=INTC NI*RNDCO)+1) 

LET L=INTCNI*RND(B)+41) 

IF B(I,J)=@ THEN 4810 

IF I#K THEN 4898 

IF J#L THEN 4698 

GOTO 4816 

LET BEI, JIJ=Bl(I,J}-! 

LET BIK,L)=BEK,L #1 

REM END OF SUBROUTINE 

RETURN 

END 
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EXERCISE 17 — (Continued) 


RUN 
INPUT SIDE OF SQUARE ARRAY?28 


INPUT ANY INTEGER BETWEEN $ AND 106720 


AFTER 8 MOVES THE RESULT IS 
QUANTA PER FREQUENCY 
SITE 


I a 


AFTER 50  #MOVES THE RESULT IS 


QUANTA PER FREQUENCY 
SITE 

8 45 

! S12 

2 41 

3 2 


AFTER 180 MOVES THE RESULT IS 


QUANTA PER FREQUENCY 
SITE 

8 77 

1 253 

2 63 

3 7 


A number of printouts have been omited. The final one is shown below. 


AFTER 508 MOVES THE RESULT IS 


QUANTA PER FREQUENCY 
SITE 

8 176 

i 168 

2 72 

3 58 

4 12 

3 2 


EXERCISE 19 — New Rules 7 


This exercise requires careful thought to decide what moves are possible 
for each oscillator as it is selected. The logic should be worked out in 
flowchart form to insure that the algorithm is correct. 





EXERCISE 21 — Goodness Of Fit 
Let f,, f2, f3,° * * be the observed frequencies. 


Let f,, f2, f3, °° * be the expected frequencies according to the Boltzmann 
distribution. 


Compute 





(Ff) 
x7 =z 


if 


then proceed as described in statistics texts. 


CHAPTER FOUR 
EXERCISE 23 — Equally-Spaced Level System 


The results are graphed below. Note that for 30 equally-spaced energy 
levels, the graph is approaching the linear relationship for an infinite 
number of energy levels. 


E/e 


(c) 
(b) 
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EXERCISE 25 — Randomly-Spaced Level System 
The results are completely determined by which energy levels happen to 
be selected. The results below reflect a program selection of the following 


levels: 2.92092, 1.35865, 4.36366, 1.45409, and 1.99161. 


E/e 


KT 
1 2 3 . 
EXERCISE 27 — Heat Capacity, Equally-Spaced Level System 
A graph of the results follows. 
Cy/k 
1.0 (c) 
(b) 
0.5 
(a) 
kT 





EXERCISE 29 — Heat Capacity, Randomly-Spaced Level System 


The same energy levels were used as in Exercise 25. A graph of the results 
follows. 


Cy/k 


0.4 


0.3 


0.2 


0.1 


EXERCISE 31 — Two Systems, Different Temperature 


System A has a higher average energy per oscillator than System B. When 
the two systems are brought into contact, quanta will be transferred back 
and forth between the two systems until the overall distribution is the 
Boltzmann distribution. The final temperature will be between Ta and 
Tp. Heat flow will take place from System A to System B. 


CHAPTER FIVE 


EXERCISE 33 — Entropy Calculation 7 


1968 REM EXERCISE 33 

{0108 REM ENTROPY SIMULATION 

18020 REM MAX SIZE YS 16 BY 18 

1030 DIM BI18,19),C(19,18),F(30} 

1935 DEF FNACX) =€X+.5)*LOGC(X) -X+.918939 
1040 PRINT “INPUT SIDE OF SQUARE ARRAY"; 
1650 INPUT Ni 

1960 MAT C=CONENI,N1) 

1865 MAT B=C€])*C 

1078 REM CARRY OUT TALLY 

1288 GOSUB 2888 

1100 REM RANDOMIZE RANDOM NUMBER GENERATOR 
111@ PRINT 

1128 PRINT “INPUT ANY INTEGER BETWEEN ! AND 168" 
1136 INPUT A 

114@ FOR I=1l TOA 

1158 LET X=RNDC@) 

1160 NEXT I 

1178 PRINT 

1175 PRINT 

1188 PRINT “MOVES,” ENTROPY/ Kk" 

1185 PRINT 

1208 REM ENTER GAME LOOPS 

1218 FOR Miz! TO 11 

1228 REM PRINT TALLY 

1230 GOSUB 386828 

1248 FOR Nz! To 25 

1250 REM MAKE A MOVE 

1268 GOSUB 49806 
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EXERCISE 33 — (Continued) 


1276 NEXT N 

1288 REM CARRY OUT TALLY 
1298 GOSUB 20082 

1388 NEXT Mi 

1318 REM END OF GAME LOOP 
1326 STOP 


2660 REM TELLY SUBROUTINE 
2018 MAT F=ZER 

2820 FOR I=! TO Ni 

2030 FOR J=1 TO Ni 

2648 LET F(B(I,J}+1 J=F(BCI,J}+1 H! 
2058 NEXT J 

2066 NEXT I 

2676 REM END OF TALLY SUBROUTINE 
286868 RETURN 

30286 REM OUTPUT SUBROUTINE 
3918 LET S=FNACNI¥NI) 

3828 FOR I=1 TO 30 

3038 j4&IF FII}3=@ THEN 5050 
3849 LET S=S-FNACF{I)) 

3058 NEXT I 

3868 PRINT CMl<1)*25,5 
3876 REM END OF SUBROUTINE 
3888 RETURN 

4000 REM SUBROUTINE FOR MOVES 
4016 LET ISI NTCNI*¥RNDCO) +1) 
4926 LET J=INTCNI*¥RNDCO) +1) 
4830 LET K=I NTCN1*RND(O)41) 
4040 LET L=INTCNI*¥RNDCO)+1) 
405@ IF BII,J)=@ THEN 4818 
4666 IF I#K THEN 48986 

40970 IF J#L THEN 4998 

4988 GOTO 4818 

4990 LET B(I,JIJ=BfI,J]~! 
4100 LET B(K,LJ=B(K,L 1 
4110 REM END OF SUBROUTINE 
4120 RETURN 

9999 END 


a) Equilibrium has been reached by about 150 moves. 


RUN 

MOVES ENTROPY/ K 
6 8 

25 935.5125 
58 129.814 
75 137.848 
188 152.544 
125 155.916 
136 162.128 
175 165.19 
2828 167.162 
225 178.5356 


258 169.695 
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EXERCISE 33 — (Continued) 


b) Equilibrium has been reached by about 75 moves. 


RUN 

MOVES ENTROPY/ K 
8 8 

25 94.5755 
28 113.294 
73 119.595 
182 1252244 
125 126.288 
158 126.299 
175 123.875 
268 127.685 
225 125.428 


258 124.522 


